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Abstract. We consider the 2D wavelet transform with 
two scales to study sky maps of temperature anisotropies 
in the cosmic microwave background radiation (CMB). 
We apply this technique to simulated maps of small sky 
patches of size 12.8° x 12.8° and 1.5' x 1.5' pixels. The 
relation to the standard approach, based on the C[s, is es- 
tablished through the introduction of the scalogram. We 
consider temperature fluctuations derived from standard, 
open and flat-A Cold Dark Matter (CDM) models. We 
analyze CMB anisotropies maps plus uncorrelated Gaus- 
sian noise (uniform and non-uniform) at different S/N lev- 
els. We explore in detail the denoising of such maps and 
compare the results with other techniques already pro- 
posed in the literature. Wavelet methods provide a good 
reconstruction of the image and power spectrum. More- 
over, they are faster than previously proposed methods. 

Key words: Cosmology: CMB - data analysis 



1. Introduction 

Future Cosmic Microwave Background (CMB) experi- 
ments will provide high resolution sky maps covering a 
wide range of frequencies. In addition to the cosmological 
CMB signal those maps will contain instrumental noise 
and contributions from Galactic and extragalactic fore- 
grounds. The denoising of these maps as well as the sep- 
aration of the different components from the CMB signal 
are the most challenging problems for CMB cosmology. 
The final goal would be to reconstruct CMB maps trying 
not to loose structural details as well as to recover the 
radiation power spectrum with the minimum error. In a 
first approach to these problems we present in this paper 
a denoising technique based on wavelets. Previously there 
have been other works based in the use of Wiener filter 



(Tegmark and Efstathiou 1996[) an d Maximum Entropy 
Methods (Hobson et al. |l99Sj , |1999| ). The use of denoising 
methods based on wavelets have certain advantages as pro- 
viding information of the contribution of different scales, 
being computationally faster (0(N)) and not requiring it- 
erative processes. The analysis of discrete 2-dimensional 
images with wavelets can be performed following different 
approaches. The two computationally faster algorithms 
are t he ones based on Multiresolution analysis (Mallat 
1989 ) and on 2D wavelet analysis (Lemarie and Meyer 
1986), using tensor products of one dimensional wavelets. 
A study of denoising of CMB maps using the former 
method has been presented in Sanz et al. ( 1999 ). This 
method is based on a single scale and three 'details' at 
each resolution level. The 2-D wavelet method used in 
this work is based on two scales, providing therefore more 
information on different resolutions (defined by the prod- 
uct of the two scales) than the Multiresolution one. More- 
over this technique is adapted to separable wavelets. On 
the other hand, an analysis of denoising using spherical 
wavelets has been recently carried out by Tenorio et al. 
( |l999| ). 

The paper is organized as follows. In §2 we present 
some basic ideas about the continuous 2D wavelet trans- 
form. We apply 2-D wavelets to the analysis of discrete 
2D CMB anisotropy maps of small sky patches in §3 and 
the conclusions are presented in §4. 



2. 2D Continuous wavelet transform 

This section is dedicated to present the continuous form 
of the 2D wavelet approach we are later using to analyse 
discrete CMB maps. As a difference with the multiresolu- 



tion approach used in Sanz et al. (1999) the 2D wavelet 
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method provides information on many more resolution el- 
ements than the former method. Moreover, this property 
is crucial for preforming an efficient linear denoising pre- 
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serving the Gaussianity of the underlying CMB field (as 
will be discussed in §3.2). 

The continuous wavelet transform of a 2D signal 
f(x\,x 2 ) is defined as 

w(R 1 ,R 2 ,bi,b 2 ) = J dxidx 2 f{xi,x 2 ) 

x^(R ll R 2 ,bi : b 2 ;xi,x 2 ) , (1) 

1 

ty(Ri,R 2 ,bi,b 2 ;x 1 ,x 2 ) 



(x\-b\ x 2 -b 2 \ 

where w(Ri, R 2 ,bi,b 2 ) is the wavelet coefficient associated 
to the scales R\ and R 2 at the point with coordinates b\ 
and b 2 . The limits in the double integral are — oo and oo 
for the two variables, ip is the wavelet 'mother' function 
that satisfies the constraints 



dxi dx 2 ip = 0, / dxi dx 2 ip = 1 



(3) 



and the 'admissibility' condition (that allows to recon- 
struct the function /), i.e. there exists the integral 



CV ee (2nY / dki dk 



hft(fci,* 2 )| 
; \hk 2 \ 



(4) 



where tp(ki, k 2 ) represents the 2D Fourier transform of ip 
and || denotes the modulus of the complex number. 

A reconstruction of the image can be achieved with 
the inversion formula 

f(xi,x 2 ) = -J- / dRl — | db 1 db 2 w(R ll R 2 ,b 1 ,b 2 ) 
W J \Ri R 2 \ 

x^/(R 1 ,R 2 ,bi,b 2 ;x 1 ,x 2 ) . (5) 

Next, let us introduce the scalogram of a 2D signal 

al(R u R 2 ) ee (w 2 (R 1 ,R 2 ,b 1 ,b 2 )) (6) 

where () means the average value calculated on the image. 

Hereinafter, we shall consider 2D wavelets that are 
separable, i.e. ip(xi,x 2 ) — ip(xi)ip(x 2 ). In this case, 

\ip(k\, k 2 )\ = \tp(ki)\ \ip(k 2 )\ . In particular, we are inter- 
ested in the Haar, the Mexican Hat and the Daubechies 
2D-transforms that can be generated in terms of the cor- 
responding ID wavelets. For the Haar case, we find 



\m\ 



- 8 . 4 (k 
nk 2 V 4 



(7) 



with an absolute maximum at k ~ 4.7 , whereas for the 
Mexican Hat 



with a single peak at k = The corresponding formu- 
lae for the Daubechies wavelets of order N can be found 
in Ogden ( [1997 ). The last wavelets form an orthonormal 
basis with compact support, increasing regularity with N 
and vanishing moments up to order TV — 1 . 

Just as an illustration we would like to present the 
scalogram for a CMB signal generated in a Standard Cold 
Dark Matter (SCDM) model. Let us assume that the im- 
age corresponds to a realization of a random field whose 
2-point correlation function is homogeneous and isotropic: 



This is equivalent to assume that the 



Fourier components f(ki,k 2 ) satisfy 
(f(k 1 ,k 2 )t(k[,k 2 )) = P(k)5(k 1 -k[)6(k 2 -k , 2 ), (9) 

where P(k), k 2 = k\ + k 2 , is the standard Fourier power 
spectrum and () means average value over realizations of 
the field (the ergodicity of the field is assumed). So, tak- 
ing average values and using equations (1) and (6), one 
obtains the variance a^ J (Ri, R 2 ) of the wavelet coefficients 
or scalogram 

al(R u R 2 )=RiR 2 J dh dk 2 P(k)\$(kiRi, k 2 R 2 )\ 2 , (10) 

For 2D white noise, i.e. P(k) — constant, one gets that 
the scalogram, cr^,, is constant at any scale. 

On the other hand, if the field / represents the tem- 
perature anisotropy of the CMB, one can obtain 




(11) 



(8) 



From the previous equation, a^jC^RxRi represents the 
power per logarithmic scale. We remark that taking into 
account the homogeneity and isotropy of the field, the 2- 
scale dependence of the scalogram is redundant in this 
case. A more appropriate treatment in this continuous ex- 
ample would be one based on isotropic wavelets defined in 
terms of a single scale. 

In the top panel of Figure |] we have represented the 
scalogram against the two scales R\ and R 2 for SCDM us- 
ing the Haar transform. The qualitative behaviour for the 
other transforms is similar. In the bottom panel we com- 
pare the scalogram along the diagonal for standard, open 
(fl = 0.3) and a flat-A (fi = 0.3, A = 0.7) CDM mod- 
els using the Haar and the Mexican Hat transforms. The 
qualitative behaviour for the two transforms is similar: 
there is a plateau for R > 1° and a maximum dependent 
on O, corresponding to the first Doppler peak. Other sec- 
ondary maxima appearing in the Figure are related to the 
secondary peaks in the standard Ci radiation power spec- 
trum (the Cg. is given by Cg ~ P(k ~ £)). Therefore, the 
position and amplitude of the maxima that appear in the 
scalogram is model dependent, being this quantity tightly 
related to the C/'s. 
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Fig. 1. Top: scalogram for a SCDM model using the Haar 
wavelet. Bottom: scalogram along the diagonal, i.e., R\ = 
R 2 , for the SCDM (solid line), open CDM (dashed line) 
and flat-A CDM (dotted line) models using the Haar (thin 
lines, bottom of the panel) and Mexican Hat (thick lines, 
top of the panel) wavelets. 

3. Denoising of CMB maps 

3.1. 2D wavelet method on a grid 

In general for a grid of 2™ x 2™ pixels, a discretization of the 
parameters of the form: R x = 2 n ~^,bi = 2 n ~^h,R 2 = 
2™ _ -' 2 ,&2 = 2 n ~i' 2 l 2 for integer- valued j and I allows to 
introduce the 2D discrete wavelet function 

^juh,h,h(^,x 2 ) = 2^+^l 2 - n ^- n x 1 -l 1 ) 

xiP(2»- n x 2 - l 2 ) , (12) 

where ji and k denote the dilation and the translation 
indexes, respectively, satisfying < ji,j 2 < n— 1,0 < h < 



2 jl - 1,0 < l 2 < 2-? 2 - I. The resolution level is defined 
by i = jl \ 32 , corresponding to R = ^R ± R 2 = 2 n ~ j . We 
also introduce a scaling function cf) that allows to define a 
complete basis to reconstruct discrete images, 

$o,o,o,o(zi,z 2 ) = 2~ n ( f ) (2- n x 1 ) ( j>(2- n x 2 ), (13) 
T *M^ X *) = 2^/ 2 ->(2-x 1 )^(2^-"x 2 - l 2 ), (14) 
rX,o,ii,o(»i.^) = V^- n ^- n Xl - h)<j>(2- n x 2 ), (15) 

We will consider orthonormal discrete bases as the 
Haar and Daubechies ones. Denoting by A any of the pre- 
vious functions, the orthonormality condition reads: 

{^h^MMi^j'^lO = ^hi'^ho'Jhl'^hl'^ ( 16 ) 

where (/, g) denotes the scalar product of two functions 
in L 2 (R 2 ). 

The wavelet coefficients are now defined by 

Wj u j 2 ,h,h=[ dx 1 dx 2 f(x 1 ,x 2 )A juj2 j 1 ,i 2 (17) 

The image is reconstructed using the following expres- 
sion: 

f(x u x 2 )= w h,hM,h A h,j2,hM( x uX2)- (18) 

jij2hh 

A representation of the wavelet coefficients can be 
done by a square that contains small squares and rect- 
angles associated to different levels of resolution. The first 
level, representing high-resolution, is ji = j 2 = 8 (i.e. 
j = 8) that contains 65536 wavelet coefficients (each one 
constructed with 2x2 = 4 pixels for the Haar trans- 
form). The second level of resolution contains two boxes: 
ji = 8,j 2 = 7 and j\ = 7,j 2 = 8 (i.e. j = 7.5) with a total 
of 2 x 32768 wavelet coefficients (each one constructed with 
2 2 x2 = 2x2 2 = 8 pixels for the Haar transform). The lev- 
els with ji = (or j 2 = 0) contain both contributions from 
wavelet-wavelet and scaling- wavelet (or wavelet-scaling). 
Finally, the lower level of resolution is j\ — j 2 = (i.e. 
j = 0) and contains four contributions: wavelet- wavelet, 
wavelet-scaling, scaling-wavelet and scaling-scaling (this 
last one is proportional to the average value of the image 
for the Haar transform). 

3.2. Reconstruction of CMB maps and radiation power 
spectra 

In the present work we have considered simulated maps 
of size 12.8° x 12.8° square degrees, pixel 1.5' x 1.5' and 
filtered with a 4.5' FWHM Gaussian beam for a standard 
CDM model (£1 = 1). We have included non-correlated 
Gaussian noise at different levels (S/N per pixel between 
0.7 and 3 at the pixel scale), considering uniform and non- 
uniform noise. This last case is introduced to account for 
the non-uniform sampling of satellite observations. As an 
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Table 1. Reconstruction errors (%) 



, ,,,,,,, 

~ ° Dn □ 




unchanged boxes 


\ 


2g 

B 

7 soft thesholded boxes 

* a 




* SURE threshold v.^ 


a ° 


removed boxes 


a 







Fig. 2. S/N ratio (open squares) in each box for a CMB 
map (SCDM model) with uniform noise at the level 
S/N — 1 using Daubechies 4. The three regions of the 
plot show the boxes that are kept unchanged (S/N > 1.5), 
removed (S/N < 0.3) or treated with a soft thresholding 
technique (boxes in between). The soft thresholds (solid 
triangles) estimated using SURE are also plotted for the 
thresholded boxes. 



extreme case we have simulated a noise map with two 
different regions, one with S/N = 3 (approximately one 
quarter of the map) and a second one with S/N = 0.7. 

The purpose of the denosing of CMB maps is to re- 
construct the original signal map as well as the radiation 
power spectrum d. 

Wavelet decompositions are performed with the pack- 
age 2D- W developed by our group. The procedure uses 
two scales, R\ — 2 n ~~ : > 1 ,R,2 = 2 11- - 72 , n — 9 in our case. 
High values of j = (ji + j<i)/1 mean high resolution, i.e., 
small scales. We distribute the coefficients Wj 1 j 2y i ly i 2 in 
boxes corresponding to a couple {31,32), having a total 
of 81 boxes. The coefficients related to the scaling func- 
tion are not included in the analysis and they arc left un- 
touched. To perform denoising, the basic operation is the 
comparison between the wavelet coefficients dispersion of 
the signal in each box with the one of the noise. The Gaus- 
sian white noise gives the same contribution in all boxes. 
Since the signal is negligible at the highest resolutions, the 
noise dispersion can be directly estimated from the map. 
Therefore, the signal dispersion can also be estimated for 
each box. 

In Figure |2[ we plot the S /N ratio (defined in terms 
of the wavelet coefficients dispersions of signal and noise) 
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for each box for a CMB simulation with S/N = 1 in real 
space. 

For the case of uniform noise, all boxes where the sig- 
nal dominates (S/N > 1.5) are kept untouched, whereas 
those with a high level of noise (S/N < 0.3) are removed. 
On the other hand, the boxes in between are treated with a 
soft thresholding technique. Given a threshold v in terms 
of the noise coefficients dispersion (a n ), the coefficients 
\w\ > va n are rescaled as w' = w ± va n (where the +,— 
signs correspond to negative and positive values of w re- 
spectively), whereas the remaining coefficients are set to 
zero. The threshold v for each box is chosen using the 
SURE method (Donoho & Johnstone jl995| ). The thresh- 
old is obtained by minimization of an unbiased estimate of 
the expected mean squared error of the estimation of the 



signal wavelet coefficients (see for instance Ogden 1997). 
In Figure ^, the thresholds obtained with the SURE tech- 
nique are plotted for a CMB map with S/N — 1. As ex- 
pected, lower S/N levels are treated with higher thresh- 
olds, i.e., more coefficients are removed. Changing the 
range of S/N where the soft technique is applied (pro- 
viding is around S/N = 1), does not appreciably change 
these results. We have used Daubechies 4 but we obtain 
little or no variations if we adopt higher order Daubechies 
wavelets. However, the Haar transform gives worse results. 
Table 1 shows the error in the map reconstructions for dif- 
ferent S/N ratios with Gaussian uniform noise. The error 
improvement achieved with the denoising technique ap- 
plied goes from factors of 3 to 5 for S/N = 3 — 0.7. The 
four top panels of Figure I show CMB maps with only 
signal (SCDM), signal plus noise with a S/N — 1, the 
reconstructed map using wavelets and the residual one. 

Regarding non-uniform (N.U.) noise, wavelet tech- 
niques allow us to treat each location in the image sep- 
arately. At each fixed location and fixed (31,32) we cal- 
culate the dispersion of the corresponding noise wavelet 
coefficient from 500 simulations of our non-uniform noise. 
Since we consider non-uniform noise that is uncorrelated, 
the average noise dispersion is the same for all the boxes, 
as in the uniform noise case. Therefore, we can get again 
the dispersion of the signal for each (31,32) pair as well as 
the S/N ratio for each coefficient. Those coefficients with 
S/N ratio in the considered range (0.3 < S/N < 1.5) are 
treated with a soft thresholding technique, whereas the 
rest are either kept or removed depending on their S/N 
ratio. Since in presence of non-uniform noise, we cannot 
use the SURE technique to estimate the optimal threshold 
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Fig. 3. Simulated map of the cosmo- 
logical signal for the SCDM model 
(top left), signal plus uniform noise 
with S/N — 1 (top right), denoised 
map using wavelets (middle left) 
and residual map obtained from the 
CMB signal map minus the denoised 
one (middle right). For comparison 
the denoised map using Wiener fil- 
ter (bottom left) is also shown to- 
gether with the corresponding resid- 
uals (bottom right). 
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Fig. 4. Power spectrum of the orig- 
inal CMB (dashed line), CMB 
plus noise (dotted line) and recon- 
structed maps using wavelets (solid 
line) for different levels of noise. Top 
left panel corresponds to S/N — 0.7, 
top-right to S/N — 1, bottom-left 
to S/N = 2 and bottom-right to 
non-uniform noise (the non-uniform 
noise map consists in two regions of 
different S/N ratio; approximately 
one quarter of the map is at the level 
of S/N=3 and the rest at S/N 
0.7; 



v (as far as we know, work is in progress to define a gen- 
eral threshold in the case of non-uniform noise, Von Sachs 



& McGibbon 1999 ), we choose for all the thresholded co- 
efficients v = 1. This threshold is defined with respect 
to the noise dispersion in each particular coefficient. We 
have chosen this value of v because it gives a good re- 
construction when comparing with the original map, but 
the results are not very sensitive to the choice of a differ- 
ent threshold in the range v — 0.8 — 2.0. In Table 1 we 
present the error of the reconstructed map in the presence 
of non-uniform noise as considered in this work. 



Regarding the power spectrum, Ce, the denoising 
method performs very well. Figures ^ and ^| show the re- 
constructed spectrum and the relative error for £ < 2000 
for different S/N ratios (the power spectrum is obtained 
in the usual way, averaging over the Fourier modes of the 
considered map at each k). The relative error is <; 20% 
for I < 1700 in all the considered cases except for the map 
with S/N = 0.7. In this last case the error increases to 
<; 30% for the same range of I's. 

Finally, we have looked for possible non-Gaussian fea- 
tures introduced by the non-linearity of the soft thresh- 
olding technique. We have compared the probability den- 
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Fig. 5. Absolute value of the rela- 
tive errors of the CMB power spec- 
trum obtained from signal plus noise 
maps (dotted line) , wavelet denoised 
maps using the SURE threshold- 
ing technique (solid line), Wiener 
denoised maps (dashed line) and 
wavelet denoised map using a lin- 
ear method, only on the top-right 
panel (dot-dashed line). Top-left 
panel corresponds to S/N — 0.7, 
top-right to S/N — 1, bottom- left to 
S/N = 2 and bottom-right to non- 
uniform noise. 



sity function of the original and the reconstructed maps 
using a Kolmogorov-Smirnov test. Both distributions are 
compatible at similar or higher levels than the original 
map and the corresponding linear reconstruction obtained 
with the Wiener Filter technique. Moreover, not signifi- 
cant change is observed in the skewness and kurtosis of 
the original and reconstructed maps. However, the appli- 
cation of soft thresholding to the wavelet coefficients at a 
certain level, which are Gaussian distributed for a Gaus- 
sian temperature field, clearly changes their distribution 
by removing all coefficients whose absolute values are be- 
low the imposed threshold and shifting the remaining ones 
by that threshold. Therefore, this technique is introducing 
a certain level of non-Gaussianity that will depend on the 
threshold imposed and that should be taken into account 
when analysing the data. If we are mainly concern about 
preserving the Gaussian character of the reconstructed im- 
age, denoising with wavelet techniques is still possible. In- 
stead of using a soft thresholding technique, we can apply 
a linear denoising method in wavelet space. In particu- 
lar, we have removed all wavelet coefficients at boxes with 
S/N < 1 and left the rest untouched. The reconstruction 
errors get only slightly worse with this simple linear tech- 
nique as shown in Table Regarding the reconstructed 
power spectrum, this is at the same level than the SURE 
reconstruction (see top-right panel of figure ||) for all the 



considered cases. It is important to remark that the linear 
denoising method based on 2D wavelets performs much 
better than the Multiresolution one due to the larger num- 
ber of boxes corresponding to the product of the 2 scales 
considered. 



3.3. Comparison with other denoising methods 

A comparis on be tween Wiener Filter (see for instance 
Press et al. 1994 ) and wavelet techniques has also been 
performed. In relation to map reconstruction the error af- 
fecting the Wiener reconstructed maps is comparable to 
that achieved with wavelet techniques in all the consid- 
ered cases. However, in order to apply Wiener filter pre- 
vious knowledge of the signal power spectrum is required. 
In a real situation this may well not be possible. The re- 
constructed and residual maps using Wiener Filter are 
shown in Figure ^. In addition, when using Wiener Filter, 
the power spectrum of the reconstructed image is clearly 
suppressed at high £'s, giving much worse results than 
the wavelet technique. For comparison, we have plotted 
in Figure |^ the absolute value of the relative error of the 
C' e s for the reconstructions with Wiener Filter. On the 
other hand, one could recover the Cg's of the original sig- 
nal by subtracting from the power spectrum of the signal 
plus noise map the estimated power spectrum of the noise, 
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which is constant at each I. However, this method gives 
in general worse results than the wavelet techniques and 
besides cannot be used to reconstruct the image. We have 
also applied a Maximum Entropy Method to the maps 
used in this work, with the definition of entropy given by 
Hobson & Lasenby ( |1998[ ). This method provides recon- 
struction errors at similar levels as the wavelet techniques. 
However, we remark that wavelet techniques are compu- 
tationally faster 0(N) and simpler to apply (not requiring 
iterative processes) than Maximum Entropy Methods. 

4. Conclusions 

We have considered the 2D wavelet transform with two 
scales to analyse CMB maps. First of all we present the 
continuous approach to the 2D wavelet. A discrete analy- 
sis is performed for a finite image. In this scaling 
function is introduced in order to define a 2D wavelet ba- 
sis. 

The 2D wavelet technique has been applied to denoise 
CMB maps for S/N ratios ranging between 0.7 and 3. 
We have also considered the case of non-uniform Gaussian 
noise. A factor between 5 and 3 is gained for the recon- 
structed images / original signal map in relation to the 
signal plus noise maps / original signal map. Regarding 
the C' e s, the relative errors are below a 20% up to I — 1700 
for all the cases with S/N > 1. A comparison with Wiener 
Filter and Maximum Entropy Method has also been per- 
formed. The later gives comparable reconstructions to the 
wavelet method, being however slower and more compli- 
cated to apply. Wiener Filter provides reconstructed maps 
with errors comparable to the wavelet technique we pro- 
pose. However, unlike the wavelet method, Wiener filter 
requires previous knowledge of the signal power spectrum. 
Moreover, the C' e s of the denoised map obtaining by ap- 
plying Wiener filter are clearly underestimated. 

Finally, we would like to remark that linear wavelet de- 
noising, which preserves Gaussianity, gives reconstruction 
errors similar to those obtained with the non-linear soft 
thresholding techniques. Moreover, linear denoising is bet- 
ter achieved by 2D wavelets than by the multiresolution 
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